% fig6.m ——— 复刻论文 Fig.6 ———
clear; clc; close all;

%% —— 模型参数 —— 
beta   = 1;      % 忆阻耦合强度 β
xi     = 0.8;    % 膜电容阻尼 ξ
mu     = 1;      % 忆阻更新速率 μ
A      = 3.4;    % 外激励幅值 A
omega  = 1;      % 外激励角频率 ω

%% —— 自适应控制参数 —— 
delta  = 0.01;   % α 增长速率 δ
lambda = 8;      % 能量阈值 λ

%% —— 仿真设置 —— 
dt     = 0.01;      % 时间步长
Tmax   = 4000;      % 总仿真时间
N      = floor(Tmax/dt);
t      = (0:N-1)*dt;  % 时间向量

%% —— 初始条件 —— 
x0     = 0.2;
y0     = 0.1;
z0     = 0.2;
alpha0 = 1.2;
S      = [x0; y0; z0];
alpha  = alpha0;

% 预分配
X      = zeros(1, N);
H      = zeros(1, N);
Alpha  = zeros(1, N);

%% —— 主循环：RK4 + 自适应 α —— 
for k = 1 : N
    % 记录
    X(k)     = S(1);
    H(k)     = S(1)^2/(2*alpha) + S(2)^2/2;
    Alpha(k) = alpha;
    
    % 当 H > λ 时，α 增长
    if H(k) > lambda
        dalpha = delta;
    else
        dalpha = 0;
    end
    
    % RK4 同步积分 [S; α]
    % —— k1 —— 
    mu_s = A*cos(omega * t(k));
    k1S = f(S, alpha, beta, xi, mu, mu_s);
    k1a = dalpha;
    % —— k2 —— 
    S2  = S + 0.5*dt*k1S;
    a2  = alpha + 0.5*dt*k1a;
    t2  = t(k) + 0.5*dt;
    mu_s2 = A*cos(omega * t2);
    k2S = f(S2, a2, beta, xi, mu, mu_s2);
    k2a = (S2(1)^2/(2*a2) + S2(2)^2/2 > lambda) * delta;
    % —— k3 —— 
    S3  = S + 0.5*dt*k2S;
    a3  = alpha + 0.5*dt*k2a;
    t3  = t(k) + 0.5*dt;
    mu_s3 = A*cos(omega * t3);
    k3S = f(S3, a3, beta, xi, mu, mu_s3);
    k3a = (S3(1)^2/(2*a3) + S3(2)^2/2 > lambda) * delta;
    % —— k4 —— 
    S4  = S +     dt*k3S;
    a4  = alpha +     dt*k3a;
    t4  = t(k) + dt;
    mu_s4 = A*cos(omega * t4);
    k4S = f(S4, a4, beta, xi, mu, mu_s4);
    k4a = (S4(1)^2/(2*a4) + S4(2)^2/2 > lambda) * delta;
    
    % 更新
    S     = S + dt*(k1S + 2*k2S + 2*k3S + k4S)/6;
    alpha = alpha + dt*(k1a + 2*k2a + 2*k3a + k4a)/6;
end

%% —— 绘图 —— 
figure;

% (a) x(τ)
subplot(3,1,1);
plot(t, X, 'k-', 'LineWidth',1);
ylabel('x','FontSize',12);
title('(a) x','FontSize',14);
xlim([0 Tmax]); grid on;

% (b) H(τ)
subplot(3,1,2);
plot(t, H, 'g-', 'LineWidth',1);
ylabel('H','FontSize',12);
title('(b) H','FontSize',14);
xlim([0 Tmax]); grid on;
text( 0.6*Tmax, lambda+0.5, sprintf('\\langleH\\rangle=%.4f', mean(H)), ...
      'FontSize',12, 'Color','k');

% (c) α(τ)
subplot(3,1,3);
plot(t, Alpha, 'r-', 'LineWidth',1);
ylabel('\alpha','FontSize',12);
xlabel('\tau','FontSize',12);
title('(c) \alpha','FontSize',14);
xlim([0 Tmax]); grid on;

% 标注论文中 α≈1.6458
[~, idx2000] = min(abs(t - 2000));
alpha2000 = Alpha(idx2000);
hold on;
plot(2000, alpha2000, 'ko','MarkerFaceColor','k');
text(2000, alpha2000+0.02, sprintf('\\alpha=%.4f', alpha2000), 'FontSize',12);

